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Abstract. The partial Schur factorization can be used to represent several eigenpairs of a matrix 
in a numerically robust way. Different adaptions of the Arnoldi method are often used to compute 
partial Schur factorizations. We propose here a technique to compute a partial Schur factorization of 
a nonlinear eigenvalue problem (NEP). The technique is inspired by the algorithm in [8], now called 
the infinite Arnoldi method. The infinite Arnoldi method is a method designed for NEPs, and can be 
interpreted as Arnoldi's method applied to a linear infinite-dimensional operator, whose reciprocal 
eigenvalues are the solutions to the NEP. As a first result we show that the invariant pairs of the 
operator are equivalent to invariant pairs of the NEP. We characterize the structure of the invariant 
pairs of the operator and show how one can carry out a modification of the infinite Arnoldi method 
by respecting the structure. This also allows us to naturally add the feature known as locking. We 
nest this algorithm with an outer iteration, where the infinite Arnoldi method for a particular type of 
structured functions is appropriately restarted. The restarting exploits the structure and is inspired 
by the well-known implicitly restarted Arnoldi method for standard eigenvalue problems. The final 
algorithm is applied to examples from a benchmark collection, showing that both processing time 
and memory consumption can be considerably reduced with the restarting technique. 

1. Introduction. The nonlinear eigenvalue problem (NEP) will in this paper 
be used to refer to the problem to find A € C C and v £ C™\{0} such that 



where M : ft — > C nx ™ is analytic in il, which is an open disc centered at the origin. 

This problem class has received a considerable amount of attention in the litera- 
ture. See, e.g., the survey papers jTHHH] and the monographs [THIS]. The results for 
( |1.1[ ) are often (but not always) presented with some restriction of the generality of 
M, such as the theory and algorithms for polynomial eigenvalue problems (PEPs) in 
[TTJ, [J_2] UH1 H] , in particular the algorithms for quadratic eigenvalue problems (QEPs) 
[351 [HE], but also recent approaches for rational eigenvalue problems (REPs) [251 177] , 
The results we will now present are directly related to [B] where an algorithm is pre- 
sented which we here call the infinite Arnoldi method. An important aspect of the 
algorithm in this paper, and the infinite Arnoldi method, is generality. Although the 
algorithm and results of this paper are applicable to PEPs, QEPs and REPs, the 
primary goal of the paper is not to solve problems for the most common structures, 
but rather to construct an algorithm which can be applied to other, less common 
NEPs in a somewhat automatic fashion. Some less common NEPs are given in the 
problem collection [2]; there exists NEPs with exponential terms [7] and implicitly 
stated NEPs such as [2"T| . 

In this paper we will present a procedure to compute a partial Schur factorization 
in the sense of the concepts of partial Schur factorizations and invariant pairs for 
nonlinear eigenvalue problems introduced in [lOj . These concepts can be summarized 
as follows. First note that the function M is in this work assumed to be analytic 
and can always be decomposed as a sum of products of constant matrices and scalar 
nonlinearities, 



M(X)v = 0, 



(1.1) 



M(A) = Mi/ X (A) + • • • + M m / m (A), 
where fi : O — > C, i = 1, . . . , m are analytic in Q. We define 

M(y, A) := MiF/i(A) + • • • + M m Y f m (A), 



(1.2) 
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where /i(A), i = 1, . . . , m are the matrix functions corresponding to /j, which are well 
defined if cr(A) C f2. An invariant pair (Y, A) G C nxp x C pxp (in the sense of QUI 
Definition 1]) satisfies 

M(Y,A)=Q. (1.3) 

Additional appropriate orthogonality conditions for Y and A yield a consistent defini- 
tion of invariant pairs and the eigenvalues of A are solutions to the nonlinear eigenvalue 
problem. In this setting, a partial Schur factorization corresponds to a particular in- 
variant pair where A is an upper triangular matrix. 

The results of this paper are based on a reformulation of the problem of finding 
an invariant pair of the nonlinear eigenvalue problem as a corresponding problem 
formulated with a (linear) infinite dimensional operator denoted B, also used in the 
infinite Arnoldi method [S]. In [5], we presented an algorithm which can be inter- 
preted as Arnoldi's method applied to the operator B. Although the operator B maps 
functions to functions, it turns out that the algorithm can be implemented with finite- 
dimensional linear algebra operations if the Arnoldi method (for B) is started with a 
constant function. This results in a Krylov subspace consisting of polynomials. Un- 
like the polynomial setting in [5] , we will in this work consider linear combinations of 
exponentials and polynomials allowing us to carry out an efficient restarting process. 
We will show that similar to the polynomial setting [8], the Arnoldi method for B 
applied to linear combinations of polynomials and exponentials can be carried out 
with finite-dimensional linear algebra operations. 

The reformulation with the operator B allows us to adapt a procedure based on 
the Arnoldi method designed for the computation of a partial Schur factorization for 
standard eigenvalue problems. We will use a construction inspired by the implicitly 
restarted Arnoldi method (IRAM) [201 H31 EH H] ■ The construction is first outlined 
adaption is outlined in Section [3] and consists of two steps respectively given in Sec- 
tion [4] and Section [5] They correspond to carrying out the Arnoldi method for the 
operator B with a locked invariant pair and a procedure to restart it. 

We finally wish to mention that there exist restarting schemes for algorithms for 



special cases of (1.1), in particular for QEPs [28), 15]. 



2. Reformulation as infinite-dimensional operator problem. In order to 



characterize the invariant pairs of (1.1) for our setting we first need to introduce some 



notation. The function B : O -> C nx ™, will be defined by 

g (A):=M(0r lM(0) : M(A) (2.1) 
A 

for A £ £A{0} an d defined as the analytic continuation at A = 0. Note that B is also 



analytic in SI, under the condition that A = is not a solution to (1.1). We will in 



this work assume that the NEP is such that A = is not an eigenvalue. From the 



definition (2.1) we reach a transformed nonlinear eigenvalue problem 

\B(\)x = x. (2.2) 



We will also use a decomposition of B similar to the decomposition ( 1.2 ) of M. That 
is, we let 



— B 1 b 1 (X) + • ■ ■ + B m b m (X), 
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(2.3) 



where bi : Q — > <C, i = 1, . . . , m are analytic in Q. Moreover, we will use the straight- 
forward coupling of the decomposition of M, by setting 

B i =M(0)- 1 M ii bi(\) = fi{0) ~ fi{X) . (2.4) 

We will use the following notation in order to express the operator and carry out 
manipulations of the operator in a concise way. Let the differentiation operator B{-^) 
be defined by the Taylor expansion in a consistent way, i.e., 

B[ Te )ip ) {e) := mm + k\ Bl{QW{e) + ^W'W + • • • - 

where <p : C — > C" is a smooth function. We are now ready to introduce the operator 
which serves as the basis for the algorithm. 

Definition 2.1 (The operator B). Let B denote the map defined by the domain 
V{B) := {ip e C 0O (R,C n ) : E^ 5(i) (°)^ (i) (°)/( i! ) is fi mte \ and the acUon 



(Btp)(0)= / <p{O)d0 + C{<p), (2.5) 
Jo 

where 

QQ 1 / J \ 

C(^):=^-B«(0)^)(0)= B(-)<^ (0). (2.6) 



i=0 v 



Several properties of the operator B are characterized in 8\. Most importantly, 



its reciprocal eigenvalues are the solutions to (2.2) and hence also to (1.1) if A ^ 0. 
In this work we will need a more general result, characterizing the invariant pairs of 
B. 

To this end we first define the application of the operator B to block functions, 
and say that if ^ : C —> C nxp with columns given by 

*(0) = (^(0),...,^,(0)), 
then B'fy is interpreted in a block fashion, i.e., 

(B9)(9) := (BM0),---,BM8))- 

With this notation, we can now consistently define an invariant pair as a pair R) 
of the operator B, where f : C ->• C nxp and R € C pxp such that 

(B*)(0) = V(6)R. (2.7) 

The following theorem explicitly shows the structure of the function ^ and relates 



invariant pairs of the operator with invariant pairs (1.3), i.e., invariant pairs in the 
setting in [10 . 

Theorem 2.2 (Invariant pairs of B). Suppose A € <C pxp is invertible and suppose 
(vf, A -1 ) is an invariant pair of B. Then, \& can be expressed as, 

V(6) =Yexp(0A), (2.8) 

for some matrix Y € C™ xp . Moreover, given A £ C pxp and 1" G ^nxp ^/jgj-g ^\ j s 
invertible, the following statements are equivalent: 
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i) The pair (^,A ), where ^(9) := Fexp(0A), is an invariant pair of the 
operator B, i.e., 

(B9){9) = ^>(9)A~ 1 . 



ii) The pair (Y, A) is an invariant pair of the nonlinear eigenvalue problem (1.1) 
in the sense of ]1(A Definition 1], i.e., 

M(Y,A)=0. (2.9) 

Proof. By differentiating (with respect to 0) the left and right-hand side of the 
definition of an invariant pair ( |2.7[ ) and using that the action of B is integration, we 
find that ^ satisfies the matrix differential equation, 

*(0) = *'(0)A _1 . 

By multiplying by A and vectorizing the equation, we have 

(A T ®/)vec(*(0)) = -^vec(¥(0)), 
a9 

and we can form an explicit solution, 

vec(*(0)) = exp(0A T <g> I)y = (exp(0A) T ® I)y Q . 



The conclusion ( |2.8[ ) follows by reversing the vectorization and setting vec(y) = j/o- 
The equivalence between statements i) and ii) follows directly from the fact that 
M(0) is invertible (since A is invertible and A = is not an eigenvalue) and the 
application of Lemma |A.1| □ 



3. Outline of the algorithm. We now know (from Theorem 2.2) that an in- 



variant pair of the NEP ( 1.1 ) is equivalent to an invariant pair of the linear operator 
B. The general idea of the procedure we will present in later sections is inspired by 
the procedures used to compute partial Schur factorizations for standard eigenvalue 
problems with the Arnoldi method [131 [UJ 123] • We will carry out a variant of the cor- 
responding algorithm for the operator B. More precisely, we will repeat the following 
two steps. 

In the first step (described in Section [4]) we compute, in a particular way, an 
orthogonal projection of the operator B onto a Krylov subspace. The projection is 
constructed such that it possesses the feature known as locking. This here means that 
given a partial Schur factorization (or an approximation of the partial Schur factor- 
ization) the projection respects the invariant subspace and returns an approximation 
containing the invariant pair (without modification) and also approximations of fur- 
ther eigenvalues. This prevents repeated convergence to the eigenvalues in the locked 
partial Schur factorization in a robust way. 

In the literature (for standard eigenvalue problems) this projection is often com- 
puted with a variation of the Arnoldi method. More precisely, we will start the Arnoldi 
algorithm with a state containing the (locked) partial Schur factorization. The re- 
sult of the infinite Arnoldi method can be expressed as what is commonly called an 
Arnoldi factorization, 



(BF k )(6) = F k+1 (6)H k , 
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(3.1) 



where H k e C (k+1}xk is a Hessenberg matrix, F k+ i : C — > C nxk is an orthogonal 
basis of the Krylov subspace and F k is the first k columns of In this paper we 

use a common notation for Hessenberg matrices; the first k rows of the matrix H k 
will be denoted H k G C . A property of the locking feature is that the Hessenberg 
matrix H_ k has the structure 

where R = (H k )i^ € C PlXpi is an upper triangular matrix. The upper left block of 



the H_ k is called the locked part, since the first pi columns of (3.1) is the equation for 



an invariant pair (2.7| 



In the first step we show how the Arnoldi method with locking can be carried 
out if we represent the functions in the algorithm (and in the factorization ( |3.1[ )) in a 
structured way. Unlike the infinite Arnoldi method in [8] we will need to work with 
functions which are linear combinations of exponentials and polynomials. It turns out 
that, similar to [5], the action of the operator as well as the entire Arnoldi algorithm 
can be carried out with finite-dimensional arithmetic, while the use of exponentials is 
benificial also for the second step. 

In the second step (described in Section™, i.e., after computing the Arnoldi 



factorization (3.1), we process the factorization such that two types of information 



can be extracted. 



We extract converged eigenvalues from the Arnoldi factorization (3.1) and 
store those in a partial Schur factorization. Due to the locking feature, the 
updated partial Schur factorization will be of the same size as the locked part 



of (3.1 1 or larger. 

• We extract a function with favorable approximation properties for those eigen- 
values of interest, which have not yet converged. 
This information is extracted in a fashion similar to implicitly restarted Arnoldi 
(IRAM) [201 EH HH HI]- However, several modifications are necessary in order to 
restart with the structured functions. 

The two steps are subsequently iterated by starting the (locked version) of Arnoldi's 
method with the extracted function and with (the possibly larger) partial Schur fac- 
torization. Thus, nesting the infinite Arnoldi method with a restarting scheme which 
is expected to eventually converge to a partial Schur factorization. 

4. The infinite Arnoldi method with locked invariant pair. In the first 
step of the conceptual algorithm described in Section [3j we need to carry out an 
Arnoldi algorithm for B with the preservation feature that the given partial Schur fac- 
torization is not modified. In an infinite-dimensional setting, the adaption to achieve 
this feature with the Arnoldi method is straightforward by initiating the state of the 
Arnoldi method with the invariant pair. The procedure is given in Algorithm[T] where 
the basis of the invariant subspace associated with the partial Schur factorization is 
assumed to be orthogonal with respect to a given scalar product < ■, • >. 

4.1. Representation of structured functions. In later sections we will pro- 
vide a specialization of all the steps in the abstract algorithm above (Algorithm [I]) 
such that we can implement it in finite-dimensional arithmetic. The first step in the 
conversion of Algorithm [T] into a finite-dimensional algorithm is to select an appro- 
priate starting function and an appropriate finite-dimensional representation of the 
functions. 
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Algorithm 1 



Input: A partial Schur factorization of B represented by (^,R) and a function / : 
C — > C™ such that < /, / >= 1 and such that / is orthogonal to the columns of 
4" with respect to < •, • >. 

Output: An Arnoldi factorization of B represented by (vi > • • • > Vfcm»* ) anc ^ 

l: Set H PuPp — R 

2: Set (ip u ...,ip pi ) = ty 

3: Set ip pi+1 = f 

4: for k = pi + 1,.. . ,k max do 

5: V = Bfki 

6: for i = 1, . . . , fe do 

7: =< 1p,tfi > 

8: ip = ip — hi tk ipi 

9: end for 

10: /lfc + l,fc = \/< V 7 , ^ > 
11: <Pk+l=1p/hk+l,k 
12: end for 



In this work, we will consider functions which are sums of exponentials and poly- 
nomials with the structure 

<p(6) = Fe se c + (4.1) 

where F G C nxp , 5 G C pxp , c G C p and g : C -> C" is a vector of polynomials. 
Moreover, we let S be a block triangular matrix 

s =( s »' £)■ (4 - 2) 

and set Six = R^ 1 G C PiXpi where p; < p, where i? will later be chosen such that 
it is an approximation of the matrix in the Schur factorization. This structure has a 
number of favorable properties important for our situation. 



The action of B applied to functions of the type (4.1 ) can be carried out in an 
efficient way using only finite-dimensional operations. This stems from the 
property that the action of B corresponds to integration and the set of poly- 
nomials and exponentials under consideration are closed under integration. 



Algorithmic details will be given in Section 4.2 

This particular structure allows the storing and orthogonalization against 
an invariant subspace, which, according to Theorem |2.2[ has exponential 
structure. 

The structure provides a freedom to choose the blocks Si 2 and S22- This 
allows us to appropriately restart the algorithm. Due to the exponential 



structure illustrated in Theorem 2.2 it will turn out to be natural to impose 
an exponential structure on the Ritz functions in order to construct a function 
/ to be used in the restart. The precise choice of S12, S22 and Y will be further 
explained in Section [5j 
In practice we also need to store the structured functions in some fashion, prefer- 
ably with matrices and vectors. It is tempting to store the exponential part and 
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polynomial part of (4.1 1 separately, i.e., to store the exponential part with the vari- 



able Y, S and c and the polynomial part by coefficients in some polynomial basis, e.g., 
the coefficients yo, . . . , yjv-i in the monomial basis q(9) — yo + y%9 + ■ ■ ■ + un-i9 N _1 - 
Although such an approach is natural from a theoretical perspective, it is not ade- 
quate from a numerical perspective. This can be seen as follows. Note that the Taylor 



expansion of the structured function (4.1 1 is 



<p(6) = (Yc + y )+ [ j [ YSc + y 1 )» + •■• + 



1 



1 

m 



YS N c 



i N 



(N 
1 

(N + iy. 



iv. 



aJV-i 



YS N+1 c 



(4.3) 



A potential source of cancellation is apparent for the first N terms in (4.3) if the poly- 



nomial q(6) approximates — Y exp(9S)c. This turns out to be a situation appearing 
in practice in this algorithm, making the storing of the structured functions in this 
separated form inadequate. 

We will instead use a function representation where the coefficients in the Taylor 
expansion are not formed by sums. This can be achieved by replacing the first iV 
terms in (4.3 1 by new coefficients xq, . . . , xjv-i, i.e., 



(YS N+1 c)9 N+1 



<p{0) = x o +x 1 + --- + x^e"- 1 + ±(YS N c)8 N + — j— ^ 

(4.4) 

The structured functions (4.1 ) will be represented with the four variables Y £ C nxp , 
S € C pxp , c e C p , x G C 7 ^, where x T = (xg\ . . . , Note that this repre- 

sentation does not suffer from the potential cancellation effects present in the naive 



representation (4.3| 



Throughout this work we will need to carry out many manipulations of functions 



represented in the form (4.4) and we need a concise notation. Let exp^y denote the 



remainder part of the truncated Taylor expansion of the exponential, i.e. 

exp N (9S) := exp(0S) - / - ±S -^S N . 



This can equivalently be expressed as 

1 



exp N (6S) 



N+l qN+1 



(N + iy. 



(N + 2)\ 



+ ■ 



(4.5) 



(4.6) 



with 



exp^eS) := exp(6»5). 



With this notation, we can now concisely express ( |4.4[ ) with exp^ and Kronecker 
products, 



ip{6) =yex Pjv _ 1 (0S)c+ ((l,9,9 2 ,...,6 N - 1 )®I n )x 



(4.7) 



4.2. Action for structured functions. We have now (in Section 4.1) intro- 
duced the function structure and shown how we can represent these functions with 
matrices and vectors. An important component in Algorithm [T] is the action of B. We 
will now show how we can compute the action of B applied to a function given with 



the representation (4.7) 



Analogous to the definition of expjy , it will be convenient to introduce a notation 
for the remainder part of the nonlinear eigenvalue problem M after a Taylor expansion 
to order N. We define, 



Mates') :=M{Y,S)-M{0)Y--^M / {0)YS-^M"(0)YS 2 ^Af^'fOjFS* 

(4.8) 

or equivalently, 

M N (Y, S) = , 1 MW + V(0)YS N+1 + , 1 „ M( N+2 H0)YS N+2 + ■■■ . (4.9) 
y ' (N+ 1)1 (N + 2)1 w y ' 

Note that with this definition 

M-iCYiS) =M(Y,S). 
We are now ready to express the action of B applied to functions with the structure 



(4.7). Note that the construction of the new function tp + = Bp in the following result 
only involves standard linear algebra operations of matrices and vectors. 

Theorem 4.1 (Action for structured functions). Let S € C pxp and c G C p be 
given constants, where S is invertible. Suppose 

<p(9) = Y expjy.i^c + ((1, 9, 8 2 , . . . , 9^) ® J„) x (4.10) 

Then, 

<p+(6) := (Btp)(B) = Yexp N (9S)c+ + ((1,9,9 2 , . . . ,9 N ) ® I n ) x+ (4.11) 

where 

c+ = S- l c, (4.12) 



A 



Or 



+>i> 



,x+,n) = (x ,...,aTjv-i) 



(4.13) 



N 



x+,o = -M(py 



(y,s)c+ + J2m^(o)x+, 



(4.14) 



Proof. We show that p+ constructed by (4.121, (4.13) and (4.14) satisfies 

Btp = tp+, (4.15) 

by first showing that the derivative of the left and the derivative of the right-hand 
side of (4.15) are equal and then showing that they are also equal in one point 8 = 0. 

■ = exp N _ 1 {6S)S. (4.16) 



From the property (4.6), we have 
d 



PYn (gq\— 1 f)NcN+l , 1 /ijV+loiV+2 

-e X p N (Vb)- m b + {N+l y° * 



Moreover, the relation (4.13) implies that 



d6 



In) X A 



((!,' 



riJV-l\ 



(4.17) 



Note that B corresponds to integration and the left-hand side of ( 4.15[ ) is <p. The right- 
hand side can be differentiated using (4.16) and (4.17). We reach that the right-hand 
side of (4.15) is <p by using (4.12). 

We have shown that the derivative of the left and the derivative of the right-hand 
side of (4.15) are equal. 



We now evaluate (4.15) at 6 = 0. From the definition of B we have that that 



(B<p)(0) — (B(-jg)ip) (0), i.e., we wish to show that 



08*0(0) = I ) (°) - <P+<P) = x o, 



(4.18) 



when N > 0. (The relation obviously holds for N = 0.) Note that by construction 
Lp + is a primitive function of ip. From the relations between /j, bi, Mi, Bi, in (2.4) it 
follows that (4.18) is equivalent to 



«'= ( (A/^i-iAW^) j s j (U) =(3/(_ )t ;__,(() ) 



(19 



'<16< 



(4.19) 



We now consider the terms of tp + in (4.11) separately. Note that for any analytic 
function g : fl — > C, we have 

5 (^)ex Pjv (0S)) (0)=g N (S), 

where is the remainder term in the truncated Taylor expansion, analogous to expjy. 
It follows that, 



For the polynomial part of tp + we have 



11 A ; 



(Y,S)c+. (4.20) 



{mQ) ((i, e, e\ . . . , e N ) ® /„) x+ )(o) = J2 m^(o) x+ , 



(4.21) 



i=0 



Note that (M(^ )y> + )(0) is the sum of ( |4.20[ ). Hence, we have sho wn (|4.19| ) (and 

hence also ( |4.18 l) by using (4.20), ( 4.21[ ) and the definition of x +i0 in ( |4.14| . 

□ 

4.3. Scalar product and finite-dimensional specialization of Algorithm[TJ 

Since the goal is to completely specify all operations in Algorithm [l] in a finite- 
dimensional setting, we also need to provide a scalar product. In [5] we worked 
with polynomials and we defined the scalar product via the Euclidean scalar product 
on monomial or Chebyshev coefficients. The structured functions described in Sec- 
tion 4.1 are not polynomials. We can however still define the scalar product consistent 
with [8 . In this work we restrict the presentation to the consistent extension of the 
definition of the scalar products via the monomial coefficients. Given two functions 



3=0 j=0 
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we define 



< <p,ip >: 



OO 

E 

i=0 



Za Xi 



(4.22) 



It is straightforward to show that (4.22) satisfies the properties of a scalar product 



and that the sum in (4.22) is always finite for functions of the considered structure. 



The computational details for the scalar product and the orthogonalization process 



are postponed until the next section (Section 4.4). 

The combination of the above results, i.e., the choice of the representation of 



the function structure (Section 4.1), the operator action (Section 4.2) and the scalar 
product (4.22), forms a complete specialization of all the operations in Algorithm [IJ 
For reasons of numerical efficiency, we will slightly modify the direct implementation 
of the operations. 

Instead of representing the individual functions (p%, . . . , (p k of the basis (tpi, . . . , ip k ) 
we will use a block representation and denote 



F k {6) = ( Vl 



,Pk) 



Now note that variables Y and S in the function structure (4.7) are not modified in 
Theorem |4.1| and obviously not modified when forming linear combinations. Hence, 
the variables Y and S can be kept constant throughout the algorithm. This allows 



us to also use the structured representation (4.7) directly for the block function F k 
instead of individually for ipi, 
matrices C k and V k such that 



■ ,fk- In every point in the algorithm, there exist 



F k (9) = Yexp N _ 1 (0S)C k + ((1, 6, . . . , 9 N - 1 ) ® I)V k , 



(4.23) 



with an appropriate choice of N. 

The variable N defining the length of the polynomial part of the structure needs 
to be adapted during the iteration. This stems from the fact that functions <p+ and ip 
in Theorem 4.1 are represented with polynomial parts of different length (N — 1 and 
N). Hence, we need to increase N by one after each application of B. Fortunately, the 
corresponding increase of N can be easily achieved by treating the leading element 
of exponential part as an element of the polynomial part. Here, this means using the 
fact that 



F h (0)=Yexp N _ 1 (6S)C k + {(l,6,- 



■i)v k 



Yexp N (9S)C k + ((1, 9, . . . , 9 N ) ® 7) ( Y ^ c ) . (4.24) 



V k 



In this work, the starting function / will be an exponential function, and after the first 
application of B, we need to expand the polynomial part with one block consisting of 
with N — 0. Since k = pi + 1 at the first application of £>, for an iteration 



YS" C k 
N\ 



corresponding to a given fc, we need to expand the polynomial part of F k with one 



block row consisting of YS N \ Ck with N 



: k - 
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Algorithm 2 Infinite Arnoldi method with structured functions and locked pair 
[V k ,C, H k ] =inf arn_exp(c, S, Y,p h k 

max,/ 

Input: Number of iterations fc max , coefficients Y G C nxp , S G C pxp , c G C p , rep- 



rcsenting the normalized function / given by (4.25) and the locked part of the 



factorization corresponding to the invariant pair R) with ^> given by ( 4.26 1 and 
R G C PlXpi is given from the structure of S in (4.27). The functions corresponding 



to the columns of \& as well as the function / must be orthogonal. 



Output: Vfe max+ i 

(£(fcmax + l)xfc max 



t + l)nx(fc„ 



Cfc n 



PpX(fcmax + l) 



representing the factorization (4.28) 



Set H 



Set C, 



Pi,pi 
Pi+i 



= R 
= (ei 



Set VJ>,+i =empty matrix of size x (pi + 1) 

for k = pi + 1, . . . , A: max do 

Compute c+ according to (4.12) where c = 

Let x G C^-W" 1 )™ be the fcth column of Vfc 

Compute x+ t i, . . . ,x + ^- Pl -i G C n according to (4.13) 
Compute x +i o according to (4.14) with N = k — pi — 1. 



Cfe, i.e., fcth column of Cfc 



Expand Vfc with one block row: 



V, 



v k 

YS k ~ Pl ~ 



(Jls-Pl— 1)1 



10: 

11: 

12 
13 
14 



Let jT fc = 



g ^(fc+l)xfe 



[c_l ,xj_,hk ,_0\ =gr am.schmidt (c + ,x + ,C k ,V k ) 
H-k-i h k 

p 

Expand C k by setting C k +i = (C k ,c±) 
Expand V k by setting V k+ i = (V k ,x±) 
end for 



With the block structure representation (4.23) we can now specialize Algorithmnl 



for the structured functions. The finite-dimensional implementation of Algorithm yy 
is given in Algorithm [2] and visually illustrated in Figure |4.1| 

The input and output of the algorithm should be interpreted as follows. The 
variables Y, S, c specify the starting function / as well as the locked part of the 
factorization. The starting function is given by 

f(0)=Yexp(0S)c (4.25) 

and the locked part of the factorization (in Algorithm [l] denoted (^,R)) corresponds 
to 

tf(0) = ycxp^S) ( Ipi ) (4.26) 



v / 

where R G C PlXpi is defined as the inverse of the leading block of S. Recall that S is 



assumed to have the block triangular structure (4.2), i.e 



R- 1 Si 



s=[ o Z ' (427) 
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The output is a finite-dimensional representation of the factorization 

(m m J(0) = ^ max+ i(#)ff fcmax , (4.28) 
where the block function -Ffc max +i is given 

*fc m «+i(0) =Yexpk m „(OS)C kmx+1 + ((l,0,>>- ,# fe — )®In)V kmc+1 , (4.29) 
and -F/c max is the first fc max columns of -Ffc max +i- 



C = 



V = 




Step 12-13 

Figure 4.1. Visualization of the infinite Arnoldi method with structured functions (Algorithm^ 

4.4. Gram-Schmidt orthogonalization. 

4.4.1. Computing the scalar product for structured functions. For struc- 



tured functions, i.e., functions of the form (|4.4|), the consistent extension of the defi- 

(4.30) 



nition (4.22) is the following. Let 

<p(0)=Yexp N (0S)c+((l,( 



3 N\ 



In)x 



and 



Then, 



tp(6) = y cx PN (es)d + ((i, e, ■ ■ • , e N ) ® i n )z. 

d H {S i ) H Y H YS l c 



N 



<(p,i)>:=J2z?Xi+ 

i=N+l 



i=0 



(*0 



M2 



(4.31) 



(4.32) 



In practice, we can compute the scalar product by truncating the infinite sum and 
exploiting the structure of the sum. 

Lemma 4.2 (Computation of scalar product). Suppose the two functions ip : C — > 



and ip : C — > C" are given by (4.30) and ( |4.31 1. Then 

N 



<<P,l/} >= ^ zf^i + d H VLAr + i iA r max C + £jv max , 



(4.33) 



i=0 
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with 



M 



i=N 



{S l ) H Y H YS i 



provides an approximation to accuracy 



(4.34) 



e 2||S|| 2 ||e||2(iW+l) 



(4.35) 



Proof. By comparing the infinite sum (4.32) with (4.33), we can solve for £jy n 
and bound the modulus, 



\e Nm J<\\dh\\Y H Y\\ 2 \\c\\ 2 J2 ^ 



2i / oo 

^<l|d||2||r ff y|| 2 || c || 2 J2 



\s\\ 



\!=W m „ + l 



The sum in the right-hand side can be interpreted as the remainder term in the Taylor 



approximation of exp(||5||2). The bound (4.35) follows by applying Taylor's theorem. 
□ 

The lemma above has some properties important from a computational perspec- 
tive: 



The sum in ( 4.34 ) involves only matrices of size px p, i.e., it does not involve 
very large matrices, under the condition that Y H Y is precomputed. 
The matrix Wn.m defined by ( 4.34[ ) is constant if S and Y are constant. 
Hence, in combination with Algorithm [2] it only needs to be computed once 
in order to construct the Arnoldi factorization. 

An appropriate value of iV max such that £j\r malc is smaller than or comparable 
to machine precision can be computed from ||iS|| by increasing -/V max until the 
right-hand side of (4.35) is sufficiently small. 



4.4.2. Computing the Gram-Schmidt orthogonalization for structured 
functions. One step of the Gram-Schmidt orthogonalization process can be seen as 
a way of computing the orthogonal complement, followed by normalizing the result. 
When working with matrices, the process is compactly expressed as follows. Consider 
an orthogonal matrix X € C" xfe . The orthogonal complement of a vector u € C, with 
respect to the space spanned by the columns of X and the Euclidean scalar product 
is given by, 



Vh. 



where 



V 



H , 



(4.36) 



(4.37) 



In the setting of Arnoldi's method, the orthogonalization coefficients h and the norm 
of the orthogonal complement j3 needs to be returned to the Arnoldi algorithm. 



Due to the fact that the considered scalar product (4.32) is the Euclidean scalar 



product on the Taylor coefficients, we can, similar to (4.36) and (4.37), compute the 



orthogonal complement using matrices. The corresponding operations for our setting 
are presented in the following theorem. 
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Algorithm 3 Gram-Schmidt orthogonalization for the scalar product (4.32) 

[cj_, x±, h, f3] =gram_schmidt(c, x, C, V) 
Input: Vectors c G C™, x G £ y ( N + 1 ) n representing the function 

<p{0) :=Ye X p N (6S)c+((l,6,--- ,6 N )®I n )x 

and C G C nxk , V G c( Ar+1 )" xfc , representing the block function, F : C -> 

F(0)=Fex Pjv (0S)C7+((l,0,--- .O^-W 



->nxk 



whose columns are orthogonal with respect to < •,- > defined by (4.32). 
Output: Orthogonalization coefficients h G C fe ,/3 G C and vectors cj_ G C p ™ and 
x± G C < -' £+1 '' Tl representing the normalized orthogonal complement of <p, 

<p ± (6) :=Yexp N {6S)c ± + ((l,6,--- ,6 N )®I n )x ± . 



1: ft = V ff a; + ^(Wjv+i.iVnt^c), where Wti+i,N msx is given by ( |4.34 | 

2: cj_ = c - Cn 

3: = X — Vh 

4: «? = + C H (WO v+ i,i W al) 

5: if ||g|| >REORTH_TOL then 

6: Cj_ = Cj_- Cg 

7: X_L = Xj_ - Vg 

8: h = h + g 

9: end if 

10: /9 = xfxj_ + C*[(W N+ i, Nmxx C_L) 

11: C_L = Cl/£ 

12: = 



Theorem 4.3 (Orthogonal complement). Let Y G C nxp , S G C pxp , C G C pxk , 
V G C n ( w + 1 ) x ' £ oe f/ig matrices representing the block function F : C — > <C nxk , 



F(9)=Yexp N (0S)C + ((l,9,- 



) ® 



where the columns are orthonormal with respect to < •,• > defined by (4.32). Consider 
the function tp, represented by c + G C p and a; + G C n ( JV+1 ^ and defined by 



and let h G C fe . 



= Yejq> JV (0S)& f + ((l,0 ) --- ,r)®4K. 



h:=V H x + +C H \ ]T -^JV2 ) c + 



(i!) 



\i=JV+l 

Then, the function (p±, represented by the vectors 

c±=c + -Ch€ C k , x± = x + — Vh G C" (JV+1) , 

and defined by 

tp±(0) :=Yexp N (pS)c± + ((l,0,-~ ,9 N )®I n )x ± 
14 



is the orthogonal complement of (p with respect to the space span by the the columns 



of F and the scalar product < ■,■ > defined by ( 4.32 ] . 

Proof. The construction is such that (p± is a linear combination of ip and the 
columns of F (due to linearity in coefficients c and x). Remains to check that (p± is 
orthogonal to columns of F. □ 

The Gram-Schmidt process with reorthogonalization can hence be efficiently im- 
plemented with operations on matrices and vectors. This is presented in Algorithm[3j 
where we used iterative reorthogonalization [3] with (as usual) at most two steps. In 
the numerical simulations we used REORTH_TOL= ^e ma ch- 

5. Extracting and restarting. Recall the general outline described in Sec- 
tion [3] and that we have now (in the Section [4]) described the first step in detail. In 
what follows we discuss the second step. We propose a procedure to carry out some 
operations of the result of the first step, i.e., Algorithm [2j and restart it such that we 
expect that the outer iteration eventually converges to a partial Schur factorization. 

5.1. Manipulations of the Arnoldi factorization. First recall that Algo- 
rithm [5] is an Arnoldi method in a function setting and the output corresponds to an 
Arnoldi factorization, 

(BF k ){6) = F k+1 {6)H k , (5.1) 
where, the block function F k +i is given by the output of Algorithm[2]with the defintion 
F k+l {6) = Yexp k (6S)C k+1 + ((l,6,--- ,6 k ) ® I n )V k+1 , (5.2) 
and F k is the first k columns of F k +i- To ease the notation, we have denoted k = A: max . 



Although the Arnoldi factorization (5.1| is a function relation, we will now see 



that several parts of the steps for implicit restarting (cf. [5D1 H2J HH [21] ) for Arnoldi's 
method (for linear matrix eigenvalue problems) can be carried out in a similar way 
by working with functions. 

We will start by computing an ordered Schur factorization of H k 

(Rn R12 i?i3\ 
R22 R23 (5.3) 
R33J 

where R n e C p ' xp ', R 22 € C^- p ^ x(p - p ^ and R 33 e c (k -^ x ^ k ~ p ^ are upper trian- 
gular matrices. The ordering is such that the eigenvalues of Rn are very accurate 
(and from now on called the locked Ritz values), the eigenvalues of R22 are wanted 
eigenvalues (selected according to some criteria) which have not converged, and the 
eigenvalues of R33 are unwanted. 
Hence, 



Q* 



H k (Q 1 ,Q 2 ,Q 3 ) = 



( Rn R12 Ri3\ 
R22 R23 
R33 

V a l a 2 «3 J 



(5.4) 



Note that ai is a measure of the (unstructured) backward error of the corresponding 
eigenvalues of Rn and ||ai|| is often used as stopping criteria. Hence, ||ai|| will be zero 
if the eigenvalues of Rn are exact and will in general be small (or very small) relative to 
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H f: since the eigenvalues of i?n are very accurate solutions. By successive application 



of Householder reflections (see e.g. 
P2 such that 



]) we can now construct an orthogonal matrix 



pi 




p 5 



i?ll 


M 





H 




e T 




p-pi 



(5.5) 



where 



is a Hessenberg matrix. 



H := 



H 



By considering the leading two blocks and columns of (5.4) and the result of the 



Householder reflection transformation (|5.5|) we find that 
\QuQ2P2)* 



H k (Q u Q 2 P 2 ) 



Rii 


Z 





H 






a[ 


e T 
p-p 







z 

H 



■o(\\ ai \\). 



(5-6) 

These operations yield a transformation of the Arnoldi factorization where the first 
block is triangular (to order 0(||ai||)). We reach the following result, which is an 
Arnoldi factorization similar to (5.f ) but only of length p. Moreover, the Hessenberg 
matrix does not contain the unwanted eigenvalues and has a leading block which is 
almost triangular. 

Theorem 5.1. Consider an Arnoldi factorization given by (5.1| and let Fk+\{9) = 
(Fk(9), f{9)). Let Qi and Q2 represent the leading blocks in the ordered Schur decom- 
position (|5.3|) and let P2, R\\ and H_ be the result of the Householder reflections in 



(5.7) 



(5.6). Moreover, let 

G p (6) := F k (B)(Q u Q 2 P 2 ), G p+1 {8) := (G p (6), f{8)). 
Then, G p+ \ approximately satisfies the length p < k Arnoldi factorization 



(BG P ){6) = G p+1 {6) 



Rn 




Z 
H 



■ o(|K 



(5.8) 



5.2. Extraction and imposing structure. Restarting in standard IRAM for 
matrices essentially consists of assigning the algorithmic state of the Arnoldi method 
to that corresponding to the factorization in Theorem |5.1| The direct adaption of 
this procedure is not suitable in our setting due to a growth of the polynomial part of 
the structured functions. This can be seen as follows. Suppose we start Algorithm [l] 
with a constant function (as done in [8]) and carry out the construction of G p as in 
Theorem 5.1 Then, G p+ \ will be a matrix with polynomials of degree k. We hence 



need to start with a state consisting of polynomials of degree k. The degree of the 
polynomial will grow with each restart and after M restarts, the polynomials will 
be of degree Mk. The representation of this polynomial will hence quickly limit the 
efficiency of the restarting scheme. 

Instead of restarting with polynomials we will perform an explicit restart using 
Algorithm [2] with a particular choice of the input which we here denote Y, S, c. This 



choice is inspired by the factorization in Theorem |5.1 
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We will first impose exponential structure on G p in the sense that we consider a 
function G p , with the property 



and defined by 



where 



G p (0) = G p (0) 

G p (9) := G p (Q)exp(Sd) 

Rn Z 
H 



Note that we can express G p (0) explicitly from (5.7) as 



G p (0) = (V fc+1 ,iQi, V k+1A Q 2 P 2 ) =: F. 



(5.9) 



(5.10) 



(5.11) 



where Vk+i,i is the upper n X (k + l)-block of 

Assume for the moment that ||ai|| = 0. Then, the first pi columns of (5.8 1 
correspond to the definition of an invariant pair (^,R), where — G Pl (6) and 

R = Rix. From Theorem 2.2 we know that VP is of exponential structure, and imposing 
the structure as in (5.9) does not modify the function, i.e., if lldi 



0, then G Pl (9) 



G Pl (9). Hence, the first pi columns of the equation (5.8) are preserved also if we replace 
Gp{&) with G p (9). Due to the fact that ||ai|| is small (or very small) we expect that 
imposing the structure as in (5.9) gives an approximation of the pi columns of (5.8), 
i.e., 



(BG Pl )(e) ~ G Pl+1 (9) 



Ru 




(5.12) 



if 



(5.9) 



ai|| is small and equality is achieved if ||ai|| = 0. 

With the above reasoning we have a justification to use the first pi columns of 
i.e., 



fexp(^) rg« 

in the initial state for the restart. In the approximation of the p — pi last columns of 



G p by the p — pi last columns of (5.9), the Arnoldi relation in the function setting is in 
general lost, because Ritz functions only have exponential structure upon convergence. 
Therefore, we will only use the (pi + l)st column in the restart, from which the Krylov 
space will be extended again in the next inner iteration This leads us to a restart with 
the function 



which corresponds to setting c = e Pi +i and initial function 

f(0)=Yexp(§6)e pt+1 . 



By these modifications of the factorization (5.8) we have now reached a choice of 
Y given by (5.11), S given by (5.10) and c = e pi +x- This choice of variables satisfy all 
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Algorithm 4 Structured explicit restarting with locking 

[S, Y] =inf arn_re start(z , A , fc max ,p) 

Input: xo G C™, Ao representing the function 

f(6) = exp(\ 6)x , 

maximum size of subspace fc max , number of wanted eigenvalues p 
Output: S, Y such that (Y, S) represents an invariant pair 

1 



9 
10 
11 
12 
13 
14 
15 
16 



=xq, with Wq jv max is given by (4.34 1 



Normalize / by setting xq — 

with S = Xq 
Set Y = (x ,0,...,0) € C nx P. 
Set S = diag(A , 1, . . . , 1) G C pxp 
Set c = e x G C p , p ; = 
while pi < p do 

[V, C,H_ km ^] =inf arn_exp(c. S'j , A: max ) 

For every eigenvalue of iifc max classify it as, lock, wanted or unwanted, and let 
p; denote the number of locked eigenvalues. 
Compute ordered Schur factorization of 



partitioned according to (5.3) 



Compute the ai vector in (5.4) 



Compute the orthogonal matrix Pi according to (5.5) 



Compute Z and H from (5.6) 



Set Yj + i = Y and Sj + \ — S according to (5.11) and (5.10) 
Reorthogonalize the function F(9) = Yj +1 exp(9Sj + i)(ei, . 
[c, •, •, ■] =gram_schmidt(e Pi+ i, ■, C j+1 , ■) 
Set j = j + 1 
end while 



the properties necessary for the input of Algorithm [2j except the orthogonality con- 
dition. The first columns of G Pl are automatically orthogonal (at least if ||ai|| = 0). 
The (pi + l)st column will however in general not be orthogonal to G Pl , which is an 
assumption needed for Algorithm [2] It is fortunately here easily remedied by orthog- 
onalizing the function corresponding to c = e pi+ i using the function gram_schmidt, 
i.e., Algorithm [3] 

The details of this selection as well as the manipulations in Section |5.1| are sum- 
marized in the outer iteration Algorithm [4j 

Remark 5.2 (Explicit restart without locking). Note that a restart which is 
theoretically very similar to what we have here proposed, can be achieved by starting the 



infinite Arnoldi method with the function of the first column of (5.9), without taking 
the "locked part" of the factorization directly into account. Such an explicit restarting 
technique (without locking) does unfortunately have unfavorable numerical properties 
and will not be persued here. From reasoning similar to [201 we know that the first 



column of (5.9) is an approximation of an element of an invariant subspace and the 
first pi steps of Arnoldi 's method started with this vector is expected to recompute the 
pi converged Ritz vectors after pi iterations. In the (pi + l)st iteration, the Arnoldi 
vector is corrupted due to cancellation. 

6. Examples. 
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6.1. A small example of Hadeler. The nonlinear eigenvalue problem pre- 
sented in [BJ, which is available with the name hadeler in the problem collection [2], 
is given by 



M (A) = -A + (A + »fA 1 + (e A +^ - 1)A 2 



where Aj g ]R nXTl , i = 0, . . . , 2 with n = 8 and /i is a shift which we will use to select 
a point close to which we will find the eigenvalues. 

In order to apply Algorithm [2] we need to derive a formula for £+,0 in (4.141. The 
derivatives for M are straightforward to compute and we compute M^r, using (4.8 1 



and (4.9). More precisely, we use the following computational expressions, 

AiY(S + nl) 2 c+ + A 2 Y(exp(S + /jj) - I)c+, 



h(Y,S)c^ 



UN 



= -AqYc 

= AiY{S 2 + 2fiS)c + + e^A 2 Y{exp{S) - 7)c H 
= AiYiS 2 c+) + e"A 2 Y(exp(S) -I- S)c+ 

S l c + " 



A > 1 



Vi = iV+l 



In the last formula, i max is chosen such that the expression has converged to machine 
precision. Since this is not computationally expensive, we can roughly overestimate 
i max . In this example it was sufficient to take i max = 40. 

In the outer algorithm (Algorithm |4| we classified a Ritz value as converged 
(locked) when the absolute residual was smaller than 1000 x £ mac h- We selected the 
largest eigenvalues of Hk as the wanted eigenvalues. 

The convergence is illustrated for two runs in Figure [6TT] and Figure [5T2] In order to 
illustrate the similarity with implicit restarting in [23j , we also carried out the infinite 
Arnoldi method with true implicit restarting by restarting only with polynomials, 
instead of using Algorithm [4] We clearly see that at least in the beginning of the 
iteration, the convergence of Algorithm [2] is similar to the convergence of IRAM. 
Note that IRAM in this setting exhibits a growth of the basis matrix and it is hence 
considerably slower. We show the number of locked Ritz- values Table [BTT] Moreover, 
we quantify the impact of the procedure to impose the structure in the restart by 
inspecting the approximation in (5.12). We define 7 as the norm of the difference 



of the left and right-hand side of (5.12). Lemma A.l shows that this difference is 
independent of 6 and provides a computable expression. More precisely, 



7 " 



(BG Pl )(6) - G Pl+1 (9) f^ 1 



G PI (0)i?n =\\M{0)- 1 M(Y,S)S- 



(6.1) 



where we used that G pi has the structure G pi (9) = Fexp(6'i? 11 1 ). The values of 7 are 



also given in Table [6TT| They are, as expected, of the same order of magnitude as the 
locking tolerance. 

6.2. A large-scale square-root example. We considered the same example 
as in [SI Section 7.2], which is the problem called gun in the problem collection [2] 
and stems from |15j . It is currently the largest example, among those examples in the 
collection 2J which are neither polynomial eigenvalue problems nor rational eigenvalue 
problems. 
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Table 6.1 

The indicator value and number of locked Ritz values for the two runs of the example of Hadeler 
in Section \6.1\ Run 1 corresponds to Figure \6.1\ and Run 2 corresponds to Figure \6.2\ The outer 
iteration count represents the number of loops carried out in Algorithm^ 




inner iteration 

Figure 6.1. Convergence of Algorithm^ (thick) and implicitly restarted Arnoldi ^23f (thin) 
ffcmax = 20, p = 10, £t = —1) for the Hadeler example in Section \6. 1\ 



In order to focus on a particular region in the complex plane we introduce (as in 
[5]) a shift fj, and a scaling 7, for which the nonlinear eigenvalue problem is 



M (A) =A - (7 A + fi)Ai + ty / 7A + yt - a\A 2 + i^A + ~\i - a\ A 3 

where u x = and <r 2 = 108.8774 and l 2 = -1. We selected 7 = 300 2 - 200 2 and 
fi = 250 2 since this transforms the region of interest to be essentially within the unit 
circle. 

In order to compute a formula for m (4.14), wc need in particular 

M{Y, S) = A Q Y - AiYijS + fjj p )+ 

iA 2 Y^J 1 S+(p-a 2 )I p + lAzY^S+ip-aDLp 



where V Z denotes the matrix square root (principal branch). 
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IRAM 

infarn restart 




10 20 30 40 50 60 70 80 
inner iteration 



Figure 6.2. Convergence of Algorithm^ (thick) and implicitly restarted Arnoldi \23f (thin) 
(inia, = 12, p = 5, p. = 3 + 5i) for the Hadeler example in Section \6.1\ 
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Figure 6.3. Computed eigenvalues and shifts for the Hadeler example 



We will partially base the formulas on the Taylor coefficients of the square root 
Ijv (needed in the computation of xq, + in ( 4.14 1 ) . We will use 



in order to compute 



where 



a ,j 
a k ,j 



^J-i\ + Ji - a 2 - = a ,j + ctij A + a 2 j\ 2 + 



(6.2a) 
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tax — 50 


^max — 30 




nof. restarts 





l 


3 


total CPU 


35.7s 


21.0s 


23.7s 


LU decomp. 


2.1s 


2.1s 


2.1s 


gram_schmidt 


23.7s 


12.9s 


15.1s 


computing x + 


6.8s 


6.9s 


3.4s 


Memory usage 


~ 200 MB 


- 78 MB 


- 58 MB 



Table 6.2 

Consumption of memory resources and profiling times, for some choices of the restart parameter 
^max o/nd p = 10. Memory in megabytes (MB) and CPU time in seconds. 



This can be used to compute of M^v, as follows, 

o(T, S)c+ = M(Y, S)c + - M(0)Yc + (6.3a) 
i(Y,S)c+ = M (Y,S)c+- M'(0)Y(Sc+) (6.3b) 

oo 1 



11 N 

' — ' I 

i=N+l 



lA 2 lY a iA Sic + J + 1 A3 I y a ^c+ J , N>1, (6.3c) 

V i=N+l J \ i=N+l J 



Note that the sums in (6.3c) are operations with vectors of relatively small dimension 
and can be computed efficiently. We selected the number of terms i max adaptively 
such that || S ,lmax c + 1| |a imaxjfc I < e ma ch- 

This results in the following formulas which we used for the computation of x+ q 

x +fi = -Af(0) -1 (Mo(Y, S)c+), for N = 
a; +: o = -M(0) -1 (Mi(y, S)c+ + M'(0)x +A ), for N = 1, 

and for N > 1, 

(N N 
M N (Y, S)c+ + lA 2 J2 x+A a iM) + lA 3 J2 x +A a iM) 
3=1 3=1 

The matrix M(0) was factorized (with an LU-factorization) before starting the iter- 
ation, such that ^^(O)^ 1 ^ could be computed efficiently 

We first wish to illustrate that the restarting and structure exploitation can con- 
siderably reduce both memory and CPU usage. In Table [5~2| we compare runs for the 
standard version of the infinite Arnoldi method [8] (first column) with the restarting 
algorithm (Algorithm |4| for two choices of the parameter fc max . The iteration was 
terminated when p = 10 eigenvalues were found. We clearly see that for the choices of 
fcmax there is a considerable reduction in memory and some reduction in computation 
time. 

In Figure |6.4| and Figure |6.5| we illustrate that the algorithm scales reasonably 
well with p, i.e., the number of wanted eigenvalues. When we increase p, we need 
more outer iterations, but eventually the algorithm usually converges for reasonably 
large p. 
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iteration 



Figure 6.4. Convergence history for Algorithm \4\ with the example involving a square root in 
Section\6^ (p = 9) 
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iteration 



Fi gure 6.5. Convergence history for Algorithm \4\ with the example involving a square root in 
Section\6^ (p = 14J 

7. Concluding remarks. We have in this work shown how the partial Schur 
factorization of an operator B can be computed using a variation of the procedures 
to compute partial Schur factorization for matrices. Several variations of the results 
for matrices appear to be possible to adapt. Concepts like thick restarting, purging 
and other selection strategies, appear to carry over but deserve further attention. We 
also wish to point that many of the results allow to be adapted or used in other ways. 
In this paper we also presented a Taylor- like scalar product, which could, essentially 
be replaced by any suitable scalar product. 
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Appendix A. A technical lemma. 

Lemma A.l. Consider Y e C" xp and S € C pxp , where S is invertible. Let 
F{9) := Yexp(6S). Then, 

(BF)(6) - F(d)S- 1 = -M(0)- x M(y, S)S~ 1 . (A.l) 



Proof. We prove the theorem by showing that the derivative of the function 
relation (A.l) holds for any 9 and that the relation holds in one point 9 — 0. Note 
that the right-hand side of (A.l I is constant (with respect to 9) and the derivative of 
the left-hand side reduces to 

F{9) - F'(9)S- 1 = F{9) -Ycxp{9S)SS- 1 = 0, 

by definition of B and differentiation of exp(9S). 

From the definition of B and evaluation of the left-hand side of (A.l ) at 9 = we 
have, 



(BF)(0) - F^S- 1 = ( B(-)Y exp(0S) ) (0) - YS' 1 . 



(A.2) 



Note that for an analytic scalar function b 



h(-)exp(0S) ) (0) = 6(5). 



Hence, 



"d9' 



B(-)YeMOS) (0) = B 1 y(6 1 (-)exp(05))(O) 



d9' 



-B m Y(b m ( — )exp(9S))(0) 



B^YhiS) 



B m Yb m (S). 



Moreover, by using the relation between bi and fi and Mi and Bi given by (2.4) we 
have, 

B^hiS) + ■ ■ ■ + B m Yb m (S) = 

M(O)" 1 [M x Y{hW ~ MS^S- 1 + ■■■ + M m Y(f m (0)I - / m (5))^ 1 ] = 

M{0y 1 [M{Y,0)S- 1 -M{Y,S)S~ 1 ] = YS^ 1 — M(0) _1 M(Y, S)S~ 1 . (A.3) 



The proof is completed by cancelling the term YS 1 when inserting (A.3) into (A.2). 
□ 
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